Simple CJS model integrated with a growth in weight model to get phi, p, and growth estimates to compare with Neural Network CMR/growth model
In all the models below, 1 = not observed and 2 = observed in the input encounter data.
Encounter data are available here in the eh.csv file. Weight data are in weight.csv
19.0.1 modelNum 1: g(i,t) model
Growth rate (gr) model structure:
gr(t,i) <- grInt(t,i)
Model code is in ./models/cmrNN_OB/tt_growth/modelCMR_tt_growth_NN_OB_functionsToSource.R
Model is run ‘by hand’ in ./models/cmrNN_OB/tt_growth/modelCMR_tt_growth_NN_OB_makeFile.R Not using this: Targets are loaded in R/modelCMR_tt_growth_NN_OB.R and saved as modelCMR_tt_growth_NN_OB_target
In the model code, a value of 1 for z or in gamma or omega signifies the individual is alive and a value of 2 signifies the individual is dead. y[i,j] is the observed data (encounter history file).
Code
out_NN_OBmod1$modelCode
{
for (i in 1:N) {
weight[i, first[i]] ~ dnorm(3.7, sd = 1.7)
weightDATA[i, first[i]] ~ dnorm(weight[i, first[i]],
sd = 0.2)
for (t in first[i]:(last[i] - 1)) {
weight[i, t + 1] <- weight[i, t] + gr[i, t] * sampleIntervalDATA[i,
t]
weightDATA[i, t + 1] ~ dnorm(weight[i, t + 1], sd = 0.2)
}
}
for (i in 1:N) {
for (t in first[i]:(last[i] - 1)) {
gr[i, t] ~ dnorm(grIntT[t], sd = 0.33)
}
}
for (t in 1:(T - 1)) {
grIntT[t] ~ T(dnorm(0, sd = 1.5), -3.5, 3.5)
grIntTOut[t] <- ilogit(grIntT[t])
}
}
Black dots/line is estimated mass and orange dots are observed masses. The green line is the first observation and the red line is the last observation. The number in the upper right corner of each panel is the fish’s cohort.
19.0.2 modelNum 2: phi(i,t) + g(i,t), p(i,t) model
Model with phi and p for each age-in-samples (t = column in the encounter history file) and individual (i)
Probability of survival (phi) model structure:
phi(t,i) <- phiInt(t,i)
Probability of capture (p) model structure:
p(t,i) <- pInt(t,i)
Growth rate (gr) model structure:
gr(t,i) <- grInt(t,i)
Survival/growth rate interaction:
Additive
Model code is in ./models/cmrNN_OB/tt_growth/modelCMR_tt_growth_NN_OB_functionsToSource.R
Model is run ‘by hand’ in ./models/cmrNN_OB/tt_growth/modelCMR_tt_growth_NN_OB_makeFile.R Not using this: Targets are loaded in R/modelCMR_tt_growth_NN_OB.R and saved as modelCMR_tt_growth_NN_OB_target
In the model code, a value of 1 for z or in gamma or omega signifies the individual is alive and a value of 2 signifies the individual is dead. y[i,j] is the observed data (encounter history file).
Code
out_NN_OBmod2$modelCode
{
for (i in 1:N) {
weight[i, first[i]] ~ dnorm(3.7, sd = 1.7)
weightDATA[i, first[i]] ~ dnorm(weight[i, first[i]],
sd = 0.2)
for (t in first[i]:(last[i] - 1)) {
weight[i, t + 1] <- weight[i, t] + gr[i, t] * sampleIntervalDATA[i,
t]
weightDATA[i, t + 1] ~ dnorm(weight[i, t + 1], sd = 0.2)
}
}
for (i in 1:N) {
for (t in first[i]:(last[i] - 1)) {
gr[i, t] ~ dnorm(grIntT[t], sd = 2)
}
}
for (t in 1:(T - 1)) {
grIntT[t] ~ T(dnorm(0, sd = 1.5), -3.5, 3.5)
}
delta[1] <- 1
delta[2] <- 0
for (i in 1:N) {
for (t in first[i]:(last[i] - 1)) {
gamma[1, 1, t, i] <- phi[i, t]
gamma[1, 2, t, i] <- 1 - phi[i, t]
gamma[2, 1, t, i] <- 0
gamma[2, 2, t, i] <- 1
}
}
for (t in 1:(T - 1)) {
p[t] ~ dunif(0, 1)
omega[1, 1, t] <- 1 - p[t]
omega[1, 2, t] <- p[t]
omega[2, 1, t] <- 1
omega[2, 2, t] <- 0
}
for (i in 1:N) {
for (t in first[i]:(last[i] - 1)) {
logit(phi[i, t]) <- phiInt[i, t]
}
}
for (i in 1:N) {
for (t in first[i]:(last[i] - 1)) {
phiInt[i, t] ~ dnorm(phiIntT[t], sd = 2)
}
}
for (t in 1:(T - 1)) {
phiIntT[t] ~ T(dnorm(0, sd = 1.5), -3.5, 3.5)
}
for (i in 1:N) {
z[i, first[i]] ~ dcat(delta[1:2])
for (j in first[i]:(last[i] - 1)) {
z[i, j + 1] ~ dcat(gamma[z[i, j], 1:2, j, i])
y[i, j + 1] ~ dcat(omega[z[i, j + 1], 1:2, j])
}
}
}
Black dots/line is estimated mass and orange dots are observed masses. The green line is the first observation and the red line is the last observation. The number in the upper right corner of each panel is the fish’s cohort.
Black dots/line is estimated probability of survival and orange dots are encountered yes (y = 0.8)/no (y = 0). The green line is the first observation and the red line is the last observation. The number in the upper right corner of each panel is the fish’s cohort.
# CJS [t,t] - growth modelSimple CJS model integrated with a growth in weight model to get phi, p, and growth estimates to compare with Neural Network CMR/growth modelIn all the models below, 1 = not observed and 2 = observed in the input *encounter data*. Encounter data are available [here](https://drive.google.com/drive/folders/1UmPv49xL-mzUBndqjRISylNO3q3Zs-cm) in the `eh.csv` file. Weight data are in `weight.csv````{r}#| label: librariesModelsNimbleRiverNN_growth#| echo: falselibrary(getWBData)library(lubridate)library(kableExtra)library(GGally)library(nimble)library(nimbleEcology)library(MCMCvis)library(tidyverse)library(targets)``````{r}#| label: globalModelsNimbleRiverCohortNN_growth#| include: falseknitr::opts_chunk$set(warning =FALSE, message =FALSE)```### modelNum 1: g(i,t) modelGrowth rate (gr) model structure: gr(t,i) \<- grInt(t,i) Model code is in `./models/cmrNN_OB/tt_growth/modelCMR_tt_growth_NN_OB_functionsToSource.R`\Model is run 'by hand' in `./models/cmrNN_OB/tt_growth/modelCMR_tt_growth_NN_OB_makeFile.R`\*Not* using this: Targets are loaded in `R/modelCMR_tt_growth_NN_OB.R` and saved as `modelCMR_tt_growth_NN_OB_target`Model runs:\#### Retrieve model results```{r}#| label: OB_tt_growth_NN#| cache: falselibrary(targets)# Following https://oliviergimenez.github.io/bayesian-cr-workshop/worksheets/4_demo.html# modelNum =1# growth only#load('./models/cmrNN_OB/tt_growth/runsOut/tt_growth_NN_OB_mostRecent.RData')load(paste0('./models/cmrNN_OB/tt_growth/runsOut/tt_growth_NN_OB_model', modelNum, '_mostRecent.RData'))out_NN_OBmod1 <- d#Input dataeh <-tar_read(eh_OB_2002_2014_target)```#### Model codeIn the model code, a value of `1` for `z` or in *gamma* or *omega* signifies the individual is alive and a value of `2` signifies the individual is dead. `y[i,j]` is the observed data (encounter history file).```{r}#| label: OB_tt_growth_NN_codeout_NN_OBmod1$modelCode```#### Model statistics```{r, echo=FALSE, results='asis'}#| label: OB_tt_growth_NN_statsTablekable(as_tibble(out_NN_OBmod1$runData), caption ="Run statistics")paste0('Run time = ', round(out_NN_OBmod1$runTime, 3), ' ', attr(out_NN_OBmod1$runTime, "units"))```#### Plot traces```{r}#| label: OB_tt_growth_NN_plotTraces priors <-runif(out_NN_OBmod1$runData$nIter * out_NN_OBmod1$runData$nChains, 0, 1)MCMCtrace(object = out_NN_OBmod1$mcmc,#ISB = FALSE,#exact = TRUE, params =c("grIntTOut"),pdf =FALSE, priors = priors)```#### Summary data```{r}#| label: OB_tt_growth_NN_plot#| MCMCplot(object = out_NN_OBmod1$mcmc, params =c("grIntTOut"))```Summary values ```{r}#| label: OB_tt_growth_NN_SummaryZ#| (summary_tt_growth <-MCMCsummary(object = out_NN_OBmod1$mcmc, params =c("grIntTOut"), round =3) %>%rownames_to_column(var ="var"))# To get the mMCMCSummaryRMNA funcion which I adapted to deal with NAssource('./models/cmrNN_OB/tt_growth/modelCMR_tt_growth_NN_OB_functionsToSource.R') summary_tt_growth_z0 <-MCMCSummaryRMNA(object = out_NN_OBmod1$mcmc, params =c("weight"), round =3) %>%rownames_to_column(var ="var") |>mutate(ind0 =str_match(var, "([0-9]+), ([0-9]+)")[,2],t0 =str_match(var, "([0-9]+), ([0-9]+)")[,3],ind =as.numeric(ind0),t =as.numeric(t0) ) |> dplyr::select(-c(t0, ind0))```Add other variables to summary values ```{r}#| label: OB_tt_growth_NN_Summary2#| ehLong <- eh$eh |>as.data.frame() |>rownames_to_column("ind0") |>pivot_longer(starts_with("ais")) |>mutate(t =as.numeric(str_match(name, "([0-9]+)")[,1]),enc = value,ind =as.numeric(ind0) ) |> dplyr::select(-c(name, value, ind0)) firstLong <- eh$first |>as_tibble() |>rownames_to_column("ind") |>mutate(ind =as.numeric(ind)) |>rename(first = value) lastLong <- eh$last |>as_tibble() |>rownames_to_column("ind") |>mutate(ind =as.numeric(ind)) |>rename(last = value) cohortLong <- eh$cohort |>as_tibble() |>rownames_to_column("ind") |>mutate(ind =as.numeric(ind)) weightLong <- eh$weight |>as.data.frame() |>rownames_to_column("ind0") |>pivot_longer(starts_with("ais")) |>mutate(t =as.numeric(str_match(name, "([0-9]+)")[,1]),weight = value,ind =as.numeric(ind0) ) |> dplyr::select(-c(name, value, ind0)) summary_tt_growth_z <- summary_tt_growth_z0 |>left_join(ehLong) |>mutate(# meanM1 = mean - 1,# pSurv = ifelse(meanM1 == 0, 1, 1 - meanM1),enc8 =ifelse(enc ==1, 0.8, 0) ) |>left_join(firstLong) |>left_join(lastLong) |>left_join(cohortLong) |>left_join(weightLong) |>mutate(resid = weight - mean ) ``````{r}#| label: OB_tt_growth_NN_getResids resids <- summary_tt_growth_z |>group_by(ind) |>summarize(meanResid =mean(abs(resid), na.rm =TRUE),n =n() ) |>ungroup() summary_tt_growth_zR <- summary_tt_growth_z |>left_join(resids)``````{r}#| cache: FALSEojs_define(numTags_OJS =dim(eh$tags)[1])ojs_define(summary_tt_growth_zR_OJS = (summary_tt_growth_zR))```#### Plot predicted and observed masses for selected individuals Select one or more individuals: ```{ojs}viewof selectInd = Inputs.select([...Array(numTags_OJS).keys()], {label:"Which fish?",value:1,step:1,multiple:true})summary_tt_growth_zR_OJS_T =transpose(summary_tt_growth_zR_OJS)summary_tt_growth_zR_OJS_T_selected = summary_tt_growth_zR_OJS_T.filter((d) => selectInd.includes(d.ind))```Black dots/line is estimated mass and orange dots are observed masses. The green line is the first observation and the red line is the last observation. The number in the upper right corner of each panel is the fish's cohort. ```{ojs}Plot.plot({width: width,height:350,inset:10,color: {scheme:"greys" },x: { label:"Age/season combination" },y: { label:"Body mass (mg)" },marks: [ Plot.frame(), Plot.dot(summary_tt_growth_zR_OJS_T_selected, {x:"t",y:"mean" }), Plot.line(summary_tt_growth_zR_OJS_T_selected, {x:"t",y:"mean" }), Plot.dot(summary_tt_growth_zR_OJS_T_selected, {x:"t",y:"weight",fill:"orange" }), Plot.ruleX(summary_tt_growth_zR_OJS_T_selected, {x:"first",y:3,"stroke":"green" }), Plot.ruleX(summary_tt_growth_zR_OJS_T_selected, {x:"last",y:3,"stroke":"red" }), Plot.text(summary_tt_growth_zR_OJS_T_selected, Plot.selectFirst({x:9,y:4,text: d =>"Cohort = "+ d.cohort }) ), Plot.text(summary_tt_growth_zR_OJS_T_selected, Plot.selectFirst({x:9,y:8,text: d =>"Residual = "+ d.meanResid }) ) ],facet: {data: summary_tt_growth_zR_OJS_T_selected,x:"ind" }})```### modelNum 2: phi(i,t) + g(i,t), p(i,t) modelModel with phi and p for each age-in-samples (t = column in the encounter history file) and individual (i) Probability of survival (phi) model structure: phi(t,i) \<- phiInt(t,i)Probability of capture (p) model structure: p(t,i) \<- pInt(t,i)Growth rate (gr) model structure: gr(t,i) \<- grInt(t,i) Survival/growth rate interaction: Additive Model code is in `./models/cmrNN_OB/tt_growth/modelCMR_tt_growth_NN_OB_functionsToSource.R`\Model is run 'by hand' in `./models/cmrNN_OB/tt_growth/modelCMR_tt_growth_NN_OB_makeFile.R`\*Not* using this: Targets are loaded in `R/modelCMR_tt_growth_NN_OB.R` and saved as `modelCMR_tt_growth_NN_OB_target`Model runs:\#### Retrieve model results```{r}#| label: OB_tt_growth_NN2#| cache: falselibrary(targets)# Following https://oliviergimenez.github.io/bayesian-cr-workshop/worksheets/4_demo.html# modelNum =2# phi + growth#load('./models/cmrNN_OB/tt_growth/runsOut/tt_growth_NN_OB_mostRecent.RData')load(paste0('./models/cmrNN_OB/tt_growth/runsOut/tt_growth_NN_OB_model', modelNum, '_mostRecent.RData'))out_NN_OBmod2 <- d#Input dataeh <-tar_read(eh_OB_2002_2014_target)```#### Model codeIn the model code, a value of `1` for `z` or in *gamma* or *omega* signifies the individual is alive and a value of `2` signifies the individual is dead. `y[i,j]` is the observed data (encounter history file).```{r}#| label: OB_tt_growth_NN_code2out_NN_OBmod2$modelCode```#### Model statistics```{r, echo=FALSE, results='asis'}#| label: OB_tt_growth_NN_statsTable2kable(as_tibble(out_NN_OBmod2$runData), caption ="Run statistics")paste0('Run time = ', round(out_NN_OBmod2$runTime, 3), ' ', attr(out_NN_OBmod2$runTime, "units"))```#### Plot traces```{r}#| label: OB_tt_growth_NN_plotTraces2 priors <-runif(out_NN_OBmod2$runData$nIter * out_NN_OBmod2$runData$nChains, 0, 1)MCMCtrace(object = out_NN_OBmod2$mcmc,#ISB = FALSE,#exact = TRUE, params =c("grIntTOut"),pdf =FALSE, priors = priors)```#### Summary data```{r}#| label: OB_tt_growth_NN_plot2#| MCMCplot(object = out_NN_OBmod2$mcmc, params =c("phiIntTOut"))MCMCplot(object = out_NN_OBmod2$mcmc, params =c("p"))MCMCplot(object = out_NN_OBmod2$mcmc, params =c("grIntTOut"))```Summary values ```{r}#| label: OB_tt_growth_NN_SummaryZ2#| (summary_tt_growth <-MCMCsummary(object = out_NN_OBmod2$mcmc, params =c("phiIntTOut", "p", "grIntTOut"), round =3) %>%rownames_to_column(var ="var"))# To get the mMCMCSummaryRMNA funcion which I adapted to deal with NAssource('./models/cmrNN_OB/tt_growth/modelCMR_tt_growth_NN_OB_functionsToSource.R') summary_tt_growth_z0 <-MCMCSummaryRMNA(object = out_NN_OBmod2$mcmc, params =c("weight"), round =3) %>%rownames_to_column(var ="var") |>mutate(ind0 =str_match(var, "([0-9]+), ([0-9]+)")[,2],t0 =str_match(var, "([0-9]+), ([0-9]+)")[,3],ind =as.numeric(ind0),t =as.numeric(t0) ) |> dplyr::select(-c(t0, ind0)) summary_tt_z0 <-MCMCSummaryRMNA(object = out_NN_OBmod2$mcmc, params =c("z"), round =3) %>%rownames_to_column(var ="var") |>mutate(ind0 =str_match(var, "([0-9]+), ([0-9]+)")[,2],t0 =str_match(var, "([0-9]+), ([0-9]+)")[,3],ind =as.numeric(ind0),t =as.numeric(t0) ) |> dplyr::select(-c(t0, ind0))```Add other variables to summary values ```{r}#| label: OB_tt_growth_NN_Summary22#| ehLong <- eh$eh |>as.data.frame() |>rownames_to_column("ind0") |>pivot_longer(starts_with("ais")) |>mutate(t =as.numeric(str_match(name, "([0-9]+)")[,1]),enc = value,ind =as.numeric(ind0) ) |> dplyr::select(-c(name, value, ind0)) firstLong <- eh$first |>as_tibble() |>rownames_to_column("ind") |>mutate(ind =as.numeric(ind)) |>rename(first = value) lastLong <- eh$last |>as_tibble() |>rownames_to_column("ind") |>mutate(ind =as.numeric(ind)) |>rename(last = value) cohortLong <- eh$cohort |>as_tibble() |>rownames_to_column("ind") |>mutate(ind =as.numeric(ind)) weightLong <- eh$weight |>as.data.frame() |>rownames_to_column("ind0") |>pivot_longer(starts_with("ais")) |>mutate(t =as.numeric(str_match(name, "([0-9]+)")[,1]),weight = value,ind =as.numeric(ind0) ) |> dplyr::select(-c(name, value, ind0)) summary_tt_growth_z <- summary_tt_growth_z0 |>left_join(ehLong) |>mutate(# meanM1 = mean - 1,# pSurv = ifelse(meanM1 == 0, 1, 1 - meanM1),enc8 =ifelse(enc ==1, 0.8, 0) ) |>left_join(firstLong) |>left_join(lastLong) |>left_join(cohortLong) |>left_join(weightLong) |>mutate(resid = weight - mean ) summary_tt_z <- summary_tt_z0 |>left_join(ehLong) |>mutate(# meanM1 = mean - 1,# pSurv = ifelse(meanM1 == 0, 1, 1 - meanM1),enc8 =ifelse(enc ==1, 0.8, 0) ) |>left_join(firstLong) |>left_join(lastLong) |>left_join(cohortLong) ``````{r}#| label: OB_tt_growth_NN_getResids2 resids <- summary_tt_growth_z |>group_by(ind) |>summarize(meanResid =mean(abs(resid), na.rm =TRUE),n =n() ) |>ungroup() summary_tt_growth_zR <- summary_tt_growth_z |>left_join(resids)``````{r}#| cache: FALSEojs_define(numTags_OJS =dim(eh$tags)[1])ojs_define(summary_tt_growth_zR_OJS = (summary_tt_growth_zR))ojs_define(summary_tt_z_OJS = (summary_tt_z))```#### Plot predicted and observed masses for selected individuals Select one or more individuals: ```{ojs}viewof selectInd = Inputs.select([...Array(numTags_OJS).keys()], {label:"Which fish?",value:1,step:1,multiple:true})summary_tt_growth_zR_OJS_T =transpose(summary_tt_growth_zR_OJS)summary_tt_growth_zR_OJS_T_selected = summary_tt_growth_zR_OJS_T.filter((d) => selectInd.includes(d.ind))summary_tt_z_OJS_T =transpose(summary_tt_z_OJS)summary_tt_z_OJS_T_selected = summary_tt_z_OJS_T.filter((d) => selectInd.includes(d.ind))```Black dots/line is estimated mass and orange dots are observed masses. The green line is the first observation and the red line is the last observation. The number in the upper right corner of each panel is the fish's cohort. ```{ojs}Plot.plot({width: width,height:350,inset:10,color: {scheme:"greys" },x: { label:"Age/season combination" },y: { label:"Body mass (mg)" },marks: [ Plot.frame(), Plot.dot(summary_tt_growth_zR_OJS_T_selected, {x:"t",y:"mean" }), Plot.line(summary_tt_growth_zR_OJS_T_selected, {x:"t",y:"mean" }), Plot.dot(summary_tt_growth_zR_OJS_T_selected, {x:"t",y:"weight",fill:"orange" }), Plot.ruleX(summary_tt_growth_zR_OJS_T_selected, {x:"first",y:3,"stroke":"green" }), Plot.ruleX(summary_tt_growth_zR_OJS_T_selected, {x:"last",y:3,"stroke":"red" }), Plot.text(summary_tt_growth_zR_OJS_T_selected, Plot.selectFirst({x:9,y:4,text: d =>"Cohort = "+ d.cohort }) ), Plot.text(summary_tt_growth_zR_OJS_T_selected, Plot.selectFirst({x:9,y:8,text: d =>"Residual = "+ d.meanResid }) ) ],facet: {data: summary_tt_growth_zR_OJS_T_selected,x:"ind" }})```#### Plot survivalBlack dots/line is estimated probability of survival and orange dots are encountered yes (y = 0.8)/no (y = 0). The green line is the first observation and the red line is the last observation. The number in the upper right corner of each panel is the fish's cohort. ```{ojs}Plot.plot({width: width,height:350,inset:10,color: {scheme:"greys" },x: { label:"Age/season combination" },y: { label:"Probability of survival" },marks: [ Plot.frame(), Plot.dot(summary_tt_z_OJS_T_selected, {x:"t",y:"pSurv" }), Plot.dot(summary_tt_z_OJS_T_selected, {x:"t",y:"enc8",fill:"orange" }), Plot.line(summary_tt_z_OJS_T_selected, {x:"t",y:"pSurv" }), Plot.ruleX(summary_tt_z_OJS_T_selected, {x:"first",y:1,"stroke":"green" }), Plot.ruleX(summary_tt_z_OJS_T_selected, {x:"last",y:1,"stroke":"red" }), Plot.text(summary_tt_z_OJS_T_selected, Plot.selectFirst({x:11,y:1,text: d => d.cohort }) ) ],facet: {data: summary_tt_z_OJS_T_selected,x:"ind" }})```#### Output model summary data for Xiaowei```{r}#| label: OB_tt_NN_xioawei#write.csv(summary_tt, './models/cmrNN_OB/tt/dataOut/xiaowei/summary_tt.csv')#write.csv(summary_tt_z, './models/cmrNN_OB/tt/dataOut/xiaowei/summary_tt_z.csv')```